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Abstract 

The Piwi-interacting RNA (piRNA) pathway is responsible for germline specification, gametogenesis, transposon silencing, and 
genome integrity. Transposable elements can disrupt genome and its functions. However, piRNA pathway evolution and its 
adaptation to transposon diversity in the teleost fish remain unknown. This article unveils evolutionary scene of piRNA pathway 
and its association with diverse transposons by systematically comparative analysis on diverse teleost fish genomes. Selective 
pressure analysis on piRNA pathway and miRNA/siRNA (microRNA/small interfering RNA) pathway genes between teleosts and 
mammals showed an accelerated evolution of piRNA pathway genes in the teleost lineages, and positive selection on functional 
PAZ (Piwi/Ago/Zwille) and Tudor domains involved in the Pi wi-piRN A/Tudor interaction, suggesting that the amino acid sub- 
stitutions are adaptive to their functions in piRNA pathway in the teleost fish species. Notably five piRNA pathway genes evolved 
faster in the swamp eel, a kind of protogynous hermaphrodite fish, than the other teleosts, indicating a differential evolution of 
piRNA pathway between the swamp eel and other gonochoristic fishes. In addition, genome-wide analysis showed higher 
diversity of transposons in the teleost fish species compared with mammals. Our results suggest that rapidly evolved piRNA 
pathway in the teleost fish is likely to be involved in the adaption to transposon diversity. 
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Introduction 

The Piwi-interacting RNA (piRNA) pathway is a small RNA 
silencing system, which functions in germline specification, 
gametogenesis, transposon silencing, genome integrity, and 
stem cell maintenance across the animal phylogeny. piRNAs 
are endogenous small noncoding RNA molecules ranging 
from 26 to 31 nucleotides in length. Animals lacking 
piRNAs exhibit activation of transposable elements (TEs) 
and defects in gametogenesis (Carmell et al. 2007). PiwiM 
(Piwi-like 1) and Piwil2 (Piwi-like 2) are core proteins inter- 
acted with piRNAs (Aravin et al. 2006; Girard et al. 2006; 
Houwing et al. 2007, 2008) and act to repress TEs both 
transcriptionally and posttranscriptionally. Other compo- 
nents involved in piRNA pathway have also been identified 
in both vertebrates and Drosophila (Ishizu et al. 2012). TDRDs 
proteins (Tudor domain-containing proteins) can associate 
with both PiwiM and Piwil2 through their N-terminal symmet- 
rical dimethyl arginines (sDMAs) in both mouse and zebrafish 



(Chuma et al. 2006; Chen et al. 2009; Reuter et al. 2009; 
Huang, Houwing, et al. 201 1). The protein arginine methyl- 
transferase 5 (Prmt5) catalyzes methylation of the arginines 
of the sDMAs (Kirino et al. 2009; Vagin et al. 2009). 
Mitochondrial protein Pld6 (Phospholipase D member 6, as 
known as MitoPLD) is also involved in primary piRNA biogen- 
esis in the mouse germline (Huang, Gao, et al. 2011; 
Watanabe et al. 201 1). Genetic analyses have revealed that 
the mouse Ddx4 (DEAD box polypeptide 4, a.k.a. Vasa, Mvh) 
is essential for the piRNA ping-pong cycling (Kuramochi- 
Miyagawa et al. 2010). Asz1 (or Gasz) (Ankyrin repeat, 
SAM, and basic leucine zipper domain containing 1) plays a 
role in retrotransposon silencing by stabilizing Piwil2 (MILI) in 
nuage (Ma et al. 2009). Both MovlOM (Moloney leukemia 
virus 10-like 1) and Fkbp6 (FK506 binding protein 6) are in- 
volved in delivering piRNAs to the mouse Piwi proteins (Frost 
et al. 201 0; Zheng et al. 201 0; Xiol et al. 201 2). Kinesin motor 
protein Kif 1 7 (Kinesin family member 17) interacts with the 
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mouse PiwiM and acts in the nucleocytoplasmic transport in 
the chromatoid bodies (Kotaja et al. 2006). In addition, Hen1 
(HEN1 methyltransferase homolog 1) is responsible for piRNA 
stability, through methylation of piRNA 3'-end of both mouse 
and zebrafish (Kirino and Mourelatos 2007; Kamminga et al. 
2010). These piRNA pathway genes, which are evolu- 
tionary conserved across the animal phylogeny, act jointly 
to protect the genome from invasive TEs in the germline 
(Ishizu etal. 2012). 

A substantial amount of piRNAs is derived from TEs, which 
play a crucial role in silencing TEs (Aravin et al. 2007). TEs are 
selfish genetic elements that can mobilize within a genome, 
potentially leading to serious genome damage. piRNAs act as 
guardians of the genome and protect it from invasive TEs in 
the germline cell (Kalmykova et al. 2005). In the process of 
gametogenesis, derepression of even a single active TE may 
cause infertility, a disaster from a Darwinian point of view. 

The types of transposons vary dramatically in vertebrate 
species, ranging from a few to over a hundred (Malone and 
Hannon 2009). In Drosophila, there are approximately 150 
different TE types, which have different expression, replica- 
tion, and mobilization strategies. There are much more TE 
types in teleost fishes than that in mammals, whereas some 
types are absent from the mammalian genomes (Aparicio 
et al. 2002; Volff et al. 2003). Around 20 clades of retrotran- 
sposons are present in pufferfish genomes, whereas only six 
clades were detected in human and mouse genomes (Volff 
et al. 2003). Notably, fish genomes contain over 30 types of 
LINE-1 retrotransposons, each of which is active, whereas 
mammalian genomes possess only one (Furano et al. 2004). 
About 4.8% transcriptome of the platyfish is derived from TE 
sequences which represent 40 different families, suggesting 
that many TEs are still active (Schartl et al. 2013). DNA trans- 
posons in fish species are often active, whereas they are inac- 
tive in mammals (Bohne et al. 2012). Although TE family 
numbers are lower in mammals than in fish species, TE copy 
numbers are generally much higher in mammals than in fishes 
(Volff 2005). Therefore, TEs in fish genomes appear to have 
undergone a higher turnover than in the mammalian ge- 
nomes (Duvernell et al. 2004), which is similar to that observed 
in Drosophila (Eickbush and Furano 2002). 

In mouse, piRNAs are mainly derived from three classes of 
retrotransposons: LTR, LINE, and SINE (Aravin et al. 2007, 
2008). However, zebrafish piRNAs have been observed to 
match to both retrotransposons and DNA transposons 
(Houwing et al. 2007, 2008) and are divided into more than 
20 groups (Zhou et al. 2010). Diverse piRNAs may provide an 
adaptive defense to repress diversified transposons in teleosts. 
These observations prompt us to test whether high diversity 
and turnover rate of TEs in fish have driven the evolution of 
piRNA pathway genes, which is involved in TE silencing. 

On the basis of their variety and abundance, TEs present an 
imposing threat to the germline cells and their genomes. With 
regard to the germline cells, they have a special need to 



protect their genome, because genetic information must be 
faithfully transmitted to offspring. Thus, it is crucial to repress 
TE activity to ensure the integrity of the germline genomes. As 
a result, piRNA machinery was suggested to have undergone 
adaptive evolution, indicating that TE activity has driven the 
evolution of piRNA pathway (Malone and Hannon 2009). 
Thus it seems that as a result from a coevolutionary arms 
race, there is a process of reciprocal adaptation between 
piRNA machinery and TE activities. Several recent studies ob- 
served evidences of selection in the piRNA pathway in 
Drosophila (Vermaak et al. 2005; Castillo et al. 2011; 
Kolaczkowski et al. 201 1; Simkin et al. 2013). So far, studies 
on adaptive evolution of piRNA pathway have been largely 
confined to the Drosophila species. As such, adaptive evolu- 
tion of piRNA pathway remains unknown in vertebrates. 

To test a possible association between piRNA pathway evo- 
lution and TE diversity in the teleost fish species, we conducted 
a comparatively evolutionary analysis for piRNA pathway 
genes between the teleost fish species and mammals, and 
compared these genes with Ago genes that function in 
miRNA/siRNA pathway. The swamp eel (Monopterus albus) 
that belongs to the order Synbranchiformes is a protogynous 
fish, which can transform its sex from female via intersex into 
male (Zhou et al. 2001; Cheng et al. 2003). We tested 
whether the selective pressure of piRNA pathway differs be- 
tween swamp eel (a kind of hermaphroditic fish) and the 
other fish species. We found an accelerated evolution of 
piRNA pathway genes in the teleost lineage. Notably, we de- 
tected many positively selected sites distributed on functional 
domains of piRNA pathway genes on the ancestral branch of 
teleosts, suggesting adaptive evolution of piRNA pathway 
genes as organisms evolve. In addition, we observed that 
five piRNA pathway genes in the swamp eel evolved more 
rapidly than other fish species, and their high expression 
level in three types of gonads (ovary, ovotestis, and testis) 
during sex reversal suggests the crucial role of piRNA pathway 
genes during the sex reversal. To further test whether high 
diversity and turnover rate of TEs exist in fish, we observed 
various types of fish TEs in comparison with mammals, sug- 
gesting a link of TEs diversity to the rapid evolution of piRNA 
pathway in fish species. Thus, this study highlights that rapidly 
evolved piRNA pathway is likely to adapt to transposon diver- 
sity in the teleost fish lineages. 

Materials and Methods 

Data Preparation 

The coding sequences of 18 genes (PiwiH, Piwil2, Tdrdi, 
Tdrkh, Pld6, Ddx4, Asz1 , MovWH, Fkbp6, Kif17, Prmt5, 
Hen1, Ago1, Ago2, Ago3a, Ago3b, Ago4, and p-actin) in 
29 vertebrates (21 mammals and 8 fishes) were obtained 
from Ensembl (www.ensembl.org, last accessed May 31, 
2014) and National Center for Biotechnology Information 
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(NCBI; http://www.ncbi.nlm.nih.gov, last accessed May 31, 
2014) (supplementary table S1, Supplementary Material 
online). Five miRNA/siRNA pathway genes (/Ago 7, Ago2, 
Ago3a, Ago3b f and Ago4), belonging to Ago subfamily in 
the Argonaute family, and fi-actin were used as control 
genes. Two analyses were performed to ensure orthology 
among genes. First, we blasted these genes in GenBank 
(http://www.ncbi.nlm.nih.gov, last accessed May 31, 2014) 
to ensure the high hits. Second, ClustalW is used to align 
the orthologous genes (Chenna et al. 2003) and to check 
their conservation among vertebrate species. A few coding 
sequences in some species (supplementary table S1, 
Supplementary Material online) are not available in Ensembl 
or NCBI because of unfinished genome sequencing. Some 
partial sequences with large gaps or sequencing errors were 
discarded in this study. 

Phylogenetic Analysis 

The protein sequences of 30 vertebrates were aligned with 
ClustalW program (Chenna et al. 2003) and manual adjust- 
ments. The maximum-likelihood tree was constructed using 
PhyML 3.0 (Guindon and Gascuel 2003; Guindon et al. 201 0) 
with LG model (Le and Gascuel 2008). To root the tree, lam- 
prey PiwiM was used as an outgroup. The phylogeny was 
analyzed using 100 replicates. The neighbor joining (NJ) 
trees using bootstrap analysis with 1 ,000 replicates were con- 
structed using MEGA 5 (Tamura et al. 2011) with Poisson 
model. 

Synonymous and Nonsynonymous Distance 

The coding sequences of 18 genes were aligned separately 
with ClustalW (Chenna et al. 2003) and manual adjustments. 
The sequences of each gene from both mammals and teleosts 
were divided into two groups (one group: 21 mammals, and 
the other: 9 teleost fish species). Synonymous (c/ s ) and non- 
synonymous (c/ N ) were calculated using the modified Nei- 
Gojobori method (Zhang et al. 1998) implemented in MEGA 
5 (Tamura et al. 201 1). The co ratios (c/ N /c/ s ) were calculated to 
compare the evolutionary distance for each gene between 
two groups. 

Branch Model Test 

To test various selection pressure acting on different clades, 
we used branch model based on the co ratio (c/ N /c/ s ) in the 
codeml in PAML 4 package (Yang 2007). A species tree of 30 
vertebrates (from Ensembl and teleost phylogeny; Chen and 
Mayden 2010) was used to determine the evolutionary rela- 
tionship. We first tested the one-ratio model as the null hy- 
pothesis, in which the entire tree has one co ratio. Two-ratio 
model was then used (a background ratio co 0 and a different 
ratio cot on the foreground branch). The significance for the 
likelihood ratio tests (LRTs) was tested using % 2 analysis. 2 Aln L 



was used (twice the difference of log likelihood between two 
models) in this analysis. 

Detection of Positive Selection 

To detect positively selected sites, the branch-site model in the 
program codeml in PAML 4 was used (Yang 2007). We set the 
ancestral branch of teleosts as foreground and the other lin- 
eages as background branch. The test 2 in PAML 4 was used 
to compare two models. In model A, three co ratios 
(0 < co 0 < 1 , co-| = 1 , co 2 > 1) are set for foreground and two 
co ratios (0 < co 0 < 1 , = 1 ) for background. Null model pa- 
rameters are the same as in model A, whereas co 2 = 1 . LRTs 
were calculated to test the significance of differences between 
two models. If the test shows positive selection, the sites 
under selection were analyzed according to high Bayes em- 
pirical Bayes posterior probabilities (BEB PPs) (Yang et al. 
2005). All positively selected sites were aligned to the protein 
sequences and functional domains in the Pfam database 
(Coggill et al. 2008; Punta et al. 2012). 

Sliding Window Analysis 

To analyze distribution of all sites with the co values (c/ N /c/ s ) 
across the PiwiH sequence, we performed sliding window 
analysis of the co values of each interval estimated by 
SWAAP (Pride 2000) with window size 60 bp and step size 
9 bp. The PiwiH sequences from 20 mammals and 8 teleosts 
were used in the analysis. 

Homology Modeling 

Three-dimensional (3D) protein structures of four zebrafish 
proteins (PiwiM -PAZ domain, Piwil2-PAZ domain, Tdrd1-the 
third Tudor domain, Tdrkh-Tudor domain) were built using 
the SWISS-MODEL server (http://swissmodel.expasy.org, last 
accessed May 31, 2014) (Schwede et al. 2003) with solved 
3D homologous proteins as templates. The crystal structures 
of templates were obtained from the Protein Data Bank (PDB). 
The PDB ID (identity) are as follows: 307V (human PiwiM -PAZ) 
(Tian etal. 201 1), 307X (human Piwil2-PAZ) (Tian etal. 201 1), 
4B9W (mouse Tdrd1-the third Tudor domain) (Mathioudakis 
et al. 2012), 3FDR (human Tdrkh-Tudor) (Chen et al. 2009). 
The models were visualized by Swiss-PDB Viewer v3.7 and 
rendered with POV-Ray v3.6 (Persistence of Vision Ray 
Tracer, USA). The positively selected sites distributed on 
these domains were highlighted as red spheres to show 
their locations in the protein structures. 

Degenerate Polymerase Chain Reaction, Rapid 
Amplification of cDNA Ends, and Reverse Transcription 
Polymerase Chain Reaction 

The swamp eels (M. albus) were obtained from markets in 
Wuhan, China. Their sexes were confirmed by microscopic 
analysis of gonad sections (ovary, ovotestis, and testis). Total 
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RNAs were isolated from adult tissues using TRIZOL 
(Invitrogen, CA). SMART cDNAs were reversely transcribed 
from the RNAs using the SMART cDNA synthesis kit following 
protocol (Clontech, CA). Degenerate PCR was used to amplify 
the conserved regions of ten genes (PiwiH, Piwil2, Tdrdi, 
Tdrkh, Prmt5, Ago1, Ago2, Ago3a, Ago3b, and Ago4) using 
multiple degenerate primers (supplementary table S2, 
Supplementary Material online). PCR was performed in a 
20 uJ reaction mix containing 10mM Tris-HCI, pH 8.3, 
1.5mM MgCI 2 , 50 mM KCI, 150mM dNTP, 0.2 jaM each 
primer, and 1 unit of Taq DNA polymerase. Based on 
the conserved regions sequences, 3 r -RACE of these genes 
was performed using common primer CDSIII: S'-ATTCTAGA 
GGCCGAGGCGGCCGACATGd(T) 30 N_1N-3 / (N = A, G, C, 
or T; N_1 = A, G, or C) and gene-specific primers (supplemen- 
tary table S3, Supplementary Material online). 5 r -RACE of 
these genes were performed using common SMARTIII 



primer: S'-AAGCAGGGTTCAACGCAGAGTGGCCATTACGGC 
CGGG-3' and gene-specific primers (supplementary table 
S3, Supplementary Material online). All sequences were 
cloned and sequenced. Semiquantitative RT-PCR was per- 
formed for five Ago genes (/Ago 7, Ago2, Ago3a, Ago3b, 
and Ago4) and six piRNA pathway genes (PiwiH, Piwil2, 
Tdrdi, Tdrkh, Ddx4, and PrmtS). Hprtwas used as a control. 
The cDNAs from adult tissues were amplified using gene- 
specific PCR primers (supplementary table S4, 
Supplementary Material online). 

TE Analysis 

To analyze diversity of TEs in fish genomes, we combined 
homology-based and de novo approaches to identify TEs in 
genomes of six teleost fish species. Zebrafish (Danio rerio), 
tetraodon (Tetraodon nigroviridis), fugu {Takifugu rubripes), 
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Fig. 1. — Phylogenetic tree of PiwiM in 30 vertebrates. Phylogenetic analysis was performed with PhyML and the lamprey (an order of jawless fish, also a 
very ancient lineage of vertebrates) PiwiM was used as an outgroup. Bootstrap values (percentage) were estimated using bootstrap analysis with 100 
replicates and the values above 50% are shown next to the nodes. The number of amino acid substitutions per site was used to assess evolutionary distances. 
The distance scale (0.1) is shown. The accession numbers of all sequences in 30 vertebrates are listed in supplementary table S1, Supplementary Material 
online. 
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medaka (Oryzias latipes), and stickleback {Gasterosteus acu- 
leatus) genome data were downloaded from Ensembl release 
63. Swamp eel (M albus) genome (DDBJ/EMBIVGenBank 
under the accession AONE00000000) was sequenced in our 
laboratory recently. The homology-based approach utilized 
database Repbase (release 1 6.02) with RepeatMasker (version 
3.2.9) (Tarailo-Graovac and Chen 2009). The de novo 
approach utilized two prediction programs (RepeatModeler 



version 1.0.4 [Price et al. 2005] and LTR-FINDER version 
1.0.5) to build the de novo repeat libraries based on the 
genome sequences. The multicopy genes and contaminations 
were removed from the libraries. Then, the RepeatMasker was 
used again to find repeats in these repetitive sequence 
libraries. Finally, we combined all the results generated by 
these methods and analyzed the diversity of TEs in fish 
genomes. 




Fig. 2. — Nonsynonymous (c/ N ) and synonymous (d s ) nucleotide distance, and d^/d s ratio for 1 8 genes between teleosts and mammals. The first 1 2 genes 
(P/W//7, Piwil2, Tdrdl, Tdrkh, Pld6, Ddx4, Asz1 , MovWH, Fkbp6, Kif17, Prmt5, and Hen 7) are involved in piRNA pathways. Five Ago genes (Ago1, Ago2, 
Ago3a, Ago3b, and Ago4), which are involved in miRNA/siRNA pathways, and fi-actin were used as control genes. (A) Nonsynonymous (c/ N ) and synon- 
ymous (d s ) nucleotide distance between teleosts and mammals. Straight lines between two circles (c/ s ) or two triangles (c/ N ) show distance between teleosts 
and mammals. All values are mean ± SE. (B) d^/d s ratio between teleosts and mammals. Black and gray columns show d^/d s ratios of different genes in 
teleosts and mammals, respectively. Gene full names are as follows: PiwiH (Piwi-like 1), Piwil2 (Piwi-like 2), Tdrdl (Tudor domain containing 1), Tdrkh (Tudor 
and KH domain containing protein), Pld6 (Phospholipase D member 6), Ddx4 (DEAD box polypeptide 4), Aszl (Ankyrin repeat, SAM and basic leucine zipper 
domain containing 1), MovWH (Moloney leukemia virus 10-like 1), Fkbp6 (FK506 binding protein 6), Kifl7 (Kinesin family member 17), Prmt5 (Protein 
arginine methyltransferase 5), Hen 1 (HEN1 methyltransferase homolog 1), Ago! (Argonaute RISC catalytic component 1), Ago2 (Argonaute RISC catalytic 
component 2), Ago3a (Argonaute RISC catalytic component 3a), Ago3b (Argonaute RISC catalytic component 3b), and Ago4 (Argonaute RISC catalytic 
component 4). The d^/d s ratios of piRNA pathway genes are significantly larger than those of miRNA/siRNA pathway genes (P< 0.01, f-test). 
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Annotation of Repeat-Associated piRNAs 

The zebrafish piRNA sequences were from Ketting's study 
(Houwing et al. 2007) together with our data (Zhou et al. 
2010). Nonredundant piRNA sequences were used to search 
against TE superfamilies using Repeat-Masking (Kohany et al. 
2006) with default parameters in Repbase (Jurka 2000). 
Percentage for each TE superfamily was calculated with re- 
spect to total TE-associated piRNAs. A total of 97,666 zebra- 
fish piRNAs were analyzed using the annotated TE database in 
Repbase. 

Results 

piRNA Pathway Genes Evolve More Rapidly in Teleosts 
than in Mammals 

To gain insights into evolutionary dynamics of piRNA pathway 
in vertebrate phylogeny, we first performed phylogenetic 
analysis of PiwiM protein, which is a core component of 
piRNA machinery in vertebrates. Comparative analysis of 
amino acid substitution rates among 30 vertebrate species, 
including 21 mammals and 9 teleost fish species, showed 
longer branch length of the PiwiM in teleosts than in mam- 
malian species (similar results in other piRNA pathway proteins 
were also observed), suggesting higher substitution rate in fish 
branches (fig. 1). To further investigate accelerated evolution 
of piRNA pathway in teleosts, using the modified Nei- 
Gojobori method (Zhang et al. 1998), we compared the 
synonymous and nonsynonymous nucleotide substitution 
distances of 12 piRNA pathway genes together with 5 
miRNA/siRNA pathway genes between two groups (mammals 
vs. teleosts). The c/ N value differences between two groups 
were markedly higher in piRNA genes (except for Hen 7) 
than in Ago genes (fig. 2A). This suggested that piRNA 
pathway genes in teleosts undertook more nonsynonymous 
substitutions than in mammals. To analyze evolutionary rates 
in piRNA pathway genes between mammals and teleosts, we 
computed the d^d s values for each gene in two groups 
(fig. 2B). We observed that c/ N /c/ s ratios of piRNA pathway 
genes were higher than that of Ago genes (P< 0.01), suggest- 
ing piRNA pathway genes evolved more rapidly than miRNA/ 
siRNA pathway genes. In addition, the d^d s ratios of 12 
piRNA pathway genes except for Hen 7 in teleosts were mark- 
edly higher than that in mammals. These results suggested 
that evolutionary rates of piRNA pathway genes were higher 
in teleosts than in mammals. Ago3b, a teleost-specific dupli- 
cated gene of Ago3 (McFarlane et al. 201 1), is an exception. 
Gene duplication leads to a relaxation of selection and conse- 
quently an elevated rate of molecular evolution for the dupli- 
cated genes (Kondrashov et al. 2002). Therefore, the c/ N /c/ s 
ratio of Ago3b in teleosts is higher than in mammals. 

To confirm these results obtained from distance analysis, 
we further performed branch model analysis by the codeml 
algorithm from the PAML 4 package (Yang 2007). 



Table 1 

Branch model tests for piRNA pathway genes and Ago genes 



Gene 


2AL a 


P Value b 


ot> 0 c 


(0-] d 


Pi will 


712.356 


6.150E-157 


0.022 


0.165 


Piwil2 


48.940 


2.639E-12 


0.081 


0.137 


Tdrdl 


8.499 


0.003 


0.253 


0.300 


Tdrkh 


33.835 


5.999E-09 


0.126 


0.210 


Pld6 


52.876 


3.552E-13 


0.036 


0.131 


Ddx4 


76.679 


2.011E-18 


0.074 


0.175 


Asz1 


35.661 


2.348E-09 


0.135 


0.283 


MovWH 


60.794 


6.336E-15 


0.099 


0.179 


Fkbp6 


31.633 


1.862E-08 


0.094 


0.185 


Kif17 


28.705 


8.428E-08 


0.063 


0.105 


Prmt5 


30.149 


4.001 E-08 


0.031 


0.059 


Hen 7 


25.956 


3.492E-07 


0.360 


0.230 


Agol 


0.402 


0.526 


0.021 


0.019 


Ago2 


27.078 


1.954E-07 


0.022 


0.009 


Ago3a 


45.255 


1.729E-11 


0.025 


0.006 


Ago3b 


9.874 


0.001 


0.025 


0.038 e 


Ago4 


79.917 


3.905E-19 


0.018 


0.003 


/3-actin 


3.359 


0.066 


0.013 


0.018 



Note. — Higher co values are in bold. 

a Twice the difference of likelihood values between two nested models (one- 
ratio model vs. two-ratio model). 

h P values significant at the 1 % threshold after Bonferroni correction in bold. 
c co 0 indicates a ratio of the mammalian clade (background). 
d (o 1 indicates a different ratio of the teleost fish clade (foreground). 
e Ago3b is a fish specific copy of Ago3. 

Comparative analysis between the two-ratio model and 
one-ratio model showed that (d^ds) values for all piRNA 
pathway genes except Hen 7 in the teleost lineage were sig- 
nificantly larger than co 0 in mammals (table 1). In contrast, 
for Ago2, Ago3a, and Ago4 in the teleost lineage were sig- 
nificantly smaller than co 0 in the mammal lineages. The co of 
Ago1 and fi-actin have no significant difference between 
mammals and fishes. These results from branch model analysis 
were consistent with synonymous and nonsynonymous dis- 
tance analysis. These data suggested that piRNA pathway 
genes evolved more rapidly in teleosts than in mammals. 

Adaptive Evolution in Functional Domains of piRNA 
Pathway Genes 

Although co ratios of piRNA pathway genes were smaller than 
1 at whole gene level, it cannot be judged that positive selec- 
tion did not occur in piRNA pathway genes. To examine the 
possible positive selection occurred in piRNA pathway, we 
used the branch-site model (Yang 2007) to test whether 
there was positive selection on the ancestral branch of the 
teleost lineages. Branch-site model test showed evidence of 
positive selection for most piRNA pathway genes (table 2). 
There were more positively selected sites in Piwi proteins 
and TDRDs proteins than other piRNA pathway proteins, for 
example, 13,12, 9, and 7 positively selected sites with BEB PPs 
higher than 95% were detected in PiwiM, Piwil2, Tdrdl, and 
Tdrkh, respectively (table 2). 
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Table 2 

Detection of positive selection for piRNA pathway genes and Ago genes 



Gene 


2AL a 


P value b 


Estimates of the 
Parameters 
in the Modified Model A c 


Positively Selected Sites 0 * 


PiwiH 


12.997 


0.0003 


co 0 = 0.068; ©2 = 40.906 
P 0 = 0.861; P 1 = 0.052 
P 2a = 0.080; P 2b = 0.004 


73K, 158A, 259S, 285R, 297C, 357V, 363V, 448S, 517G, 531Y, 535Q, 555V, 684L, 
692S, 850N, 852D 


Piwil2 


17.351 


3.106E-05 


©o = 0.083; ©2=10.057 
P 0 = 0.826; P 1 = 0.076 
P 2a = 0.089; P 2b = 0.008 


305P, 340T, 345V, 386P, 553S, 556G, 602T, 647C, 666V, 672S, 692E, 7131, 735S, 
737P, 804T, 827S7839S7983Y 


Tdrdl 


28.065 


1.173E-07 


©o = 0.183; ©2 = 86.504 
P 0 = 0.588; Pt = 0.199 
P 2a = 0.157; P 2b = 0.053 


5F, 14L, 31P, 91N, 161Y, 163R, 352Q, 438N, 453R, 474K, 496E, 518P, 662E, 776S, 
787E, 790G, 791 Q, 793V, 922P, 1000P, 1011F, 1054L 


Tdrkh 


20.404 


6.269E-06 


©o = 0.116; ©2 = 60.256 
P o = 0.598; Pt = 0.146 
P 2a = 0.204; P 2b = 0.050 


69H, 123Q, 126V, 131Q, 138L, 142C, 152E, 157V, 159D, 177K, 190S, 198E, 21 3Q, 
241 Q, 251 P, 348S, 373E, 385Y, 387D, 402S, 404L, 41 8L, 453T, 480S, 505K, 51 5S, 
542G 


Pld6 


4.124 


0.042 


©o = 0.065; ©2 = 6.764 
P 0 = 0.784; P 1 = 0.044 
P 2a = 0.161; P 2b = 0.009 


59H, 68P, 82S, 94S, 120V, 146S, 166G, 178L, 179T, 196L 


Ddx4 


9.197 


0.002 


©o = 0.068; ©2 = 335.073 
P 0 = 0.766; P 1 = 0.177 
P 2a = 0.045; P 2b = 0.010 


51 S, 104F, 177R, 268S, 337Q, 357A 


Asz1 


10.732 


0.001 


©o = 0.155; ©2 = 40.485 
P 0 = 0.700; P 1 = 0.146 
P 2a = 0.126; P 2b = 0.026 


154C, 1901, 396K, 41 4Q, 41 9C 


MovWH 


17.078 


3.587E-05 


©o = 0.085; ©2=11.126 
P 0 = 0.694; P 1 = 0.130 
P 2a = 0.146; P 2b = 0.027 


302K, 31 1S, 341 P, 366H, 390P, 409V, 438A, 560I, 561 G # 737Q, 762V, 775H, 778F, 
779H, 787S, 866N, 875G, 886D, 892K, 908C, 91 2P, 91 3R, 920C, 921 Q 


Fkbp6 


6.917 


0.009 


©o = 0.109; ©2 = 46.815 
P 0 = 0.721; P 1 = 0.049 
P 2a = 0.213; P 2b = 0.014 


79T, 96A, 118Q, 158S, 164T 


Kif17 


8.569 


0.003 


©o = 0.045; ©2 = 4.658 
P 0 = 0.798; P 1 = 0.099 
P 2a = 0.090; P 2b = 0.011 


122F, 222S, 226A, 366G, 560L, 575S, 581S, 595H, 640G, 641Q, 645G, 651V, 678T, 
684V, 730E 


Prmt5 


19.085 


1.250E-05 


©o = 0.036; ©2 = 97.992 
P 0 = 0.874; P 1 = 0.030 
P 2a = 0.091; P 2b = 0.003 


5S, 51E, 95E f 110C, 124G, 125P, 135L, 187S, 199C, 212T, 216K, 269S, 407D, 463S, 
488K, 492C, 523D, 530Q7531C, 538C7^52f7^62K, 587S, 590D, 602G 


Hen1 


4.658 


0.031 


©o = 0.160; ©2 = 99.653 
P 0 = 0.602; P 1 = 0.351 
P 2a = 0.029; P 2b = 0.016 


179C 


Ago1 


6.078 


0.013 


©o = 0.017; ©2 = 9.523 
P 0 = 0.979; P 1 = 0.008 
P 2a = 0.011; P 2b = 0.0001 


131 R, 716S 


Ago2 


2E-06 


0.999 


©o = 0.012; ©2=1 
P 0 = 0.964; P 1 = 0.004 
P 2a = 0.030; P 2b = 0.0001 




Ago3a 


0 


1 


r\ r\ a n . /~ f\C\C\ 

©o = 0.014; ©2 = 6.090 
P 0 = 0.999; P 1 = 0.0003 

P 2a = 0; P 2b = 0 




Ago3b 


0 


1 


©o = 0.030; ©2=12.715 
P 0 = 0.992; P 1 = 0.007 

P 2a = 0; P 2b = 0 




Ago4 


0 


1 


©o = 0.007; ©2=10.684 
Po=1; ?1 = 0 
P 2a = 0; P 2b = 0 




p-actin 


4.324 


0.037 


©o = 0.012; ©2 = 998.998 
P 0 = 0.982; P 1 = 0.009 
P 2a = 0.008; P 2b = 0.00007 





a The ancestral branch leading to the teleosts was assigned as the foreground branch. Twice the difference between the likelihood values of the null model and the 
alternative model is indicated as 2AL The modified model A is used as the alternative model and the modified model A with co 2 fixed at 1 is the null model. 
b P values significant at the 5% threshold after Bonferroni correction in bold. 

C P 0 is the proportion of codons that have co 0 in all branches, P^ is the proportion of codons that have cdt = 1 in all branches, P 2a is the proportion of codons that have co 0 in the 
background branches but co 2 in the foreground branches, and P 2 b is the proportion of codons that have in the background branches but co 2 in the foreground branches. 

d Sites with the BEB PPs higher than 90% are shown. Sites higher than 95% are single underlined and sites higher than 99% are double underlined. The sites of MovlOM 
and Fkbp6 are indexed by the amino acids at the sites in the tetraodon, whereas the sites of the other genes are indexed by the amino acids at the sites in the zebraf ish. 



Genome Biol. Evol. 6(6): 1393-1 407. doi:10.1093/gbe/evu105 Advance Access publication May 19, 2014 



1399 



Yi etal. 



GBE 



A 

c 10 1 
5 




0 500 1000 1500 2000 2500 b p 

PiwiM H3SM W7*W MID PIWI 



Fig. 3. — Distribution of the positively selected sites in functional domains. (A) Putatively positively selected sites with BEB PPs higher than 90% were 
shown. Functional domains in x axis: PAZ-MID (Piwi/Ago/Zwille-Middle), PIWI (P-element induced wimpy testis), Tudor, KH (K homology), PLDc 
(Phospholipase D family), DEAD (Asp-Glu-Ala-Asp box), Ank (Ankyrin repeat), DEXD (DEAD-like helicase superfamily ATP binding domain), AAA12 
(ATPases Associated with a wide variety of cellular Activities superfamily), FKBPc (FKBP-type peptidyl-prolyl cis-trans isomerase), Kinesin (Kinesin motor 
domain), BBP1c (C -terminal domain of BBP1), Mtase (Arginine-A/-methyltransferase). The numbers above stacked columns show the ratios of the sites 
distributed on functional domains/total positively selected sites for each gene. (B) Sliding window analysis shows co value distribution (c/n/c/ s ) along the coding 
sequence of PiwiH between teleosts (red line) and mammals (blue line). Functional domains of PiwiM : PAZ, MID, and PIWI. sDMA represents the binding sites 
for Tudor domain of Tdrdl or Tdrkh. 



We further analyzed locations of the positively selected 
sites in these genes. Using the protein topologies, we ob- 
served that their distributions were clearly concentrated on 
functional domains (fig. 3A). As for Piwi proteins, 50% 
(PiwiM) and 56% (Piwil2) positively selected sites were distrib- 
uted in PAZ-MID (Piwi/Ago/Zwille-Middle) domain, which is 



responsible for recognition of piRNA 3'-/5'-end (Tian et al. 
2011). In addition, 31% (PiwiM) and 22% (Piwil2) positively 
selected sites were distributed in PIWI domain, which is similar 
to ribonuclease H (Song et al. 2004; Tolia and Joshua-Tor 
2007). Branch-site results were further supported by sliding 
window analysis, which showed variation in co values along 
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Fig. 4. — Structural modeling suggests functional importance of positively selected sites. NJ trees recapitulate vertebrate phylogeny. The amino acid 
alignments of nine representative species are shown in the middle panels. Homology modeling of zebrafish PiwiM PAZ domain (from 275E to 387T) (A), 
zebrafish Piwil2 PAZ domain (from 459R to 574T) (B), the third Tudor domain (from 61 6L to 807G) of zebrafish Tdrdl (0, and zebrafish Tdrkh Tudor domain 
(from 348S to 440C) (D). The binding pockets (piRNA 3'-end binds to PAZ; sDMA of Piwi binds to Tudor) are shown with dashed ellipse. Arrows (in amino 
acid alignments) or red spheres (in 3D models) indicate the positively selected sites. The modeling domains are colored by structural motifs: a helix, blue; p 
sheet, cyan; coil, purple. The modeling structures in PiwiM , Piwil2, and Tdrkh are not a whole domain because of incompleteness of corresponding template. 



the PiwiH gene (fig. 3B). These data suggested that positive 
selection was strongest at the interface between Piwi proteins 
and target piRNA molecules, implicating their possible roles in 
modulating interaction of Piwi proteins with diverse piRNAs in 
teleosts. 



As for two TDRDs proteins, 50% and 30% positively se- 
lected sites were distributed in Tudor domain of Tdrdl and 
Tdrkh, respectively. Tudor domain can interact with sDMA in 
the N-terminal region of Piwi proteins (Chen et al. 2009; 
Mathioudakis et al. 201 2). In addition, 44% positively selected 
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Fig. 5. — Branch model analysis of the co values for 12 genes in the teleost lineage. (A) Cloning and expression levels of five Ago genes and six piRNA 
pathway genes in the swamp eel detected by RT-PCR. Hprt was used as an inner control. Ago genes are expressed ubiquitously, whereas the piRNA pathway 
genes are mainly expressed in three kinds of gonads except Tdrkh. (B) The branch model analysis (one-ratio model against two-ratio model). Black columns 
indicate co values of the genes in the swamp eel (foreground branch), and gray columns indicate co values of the other fish species (background branch). 
*LRT, P< 0.05; **p< 0.01 . Ago genes and p-actin were used as controls. 



sites were distributed in Tdrkh-KH domain, which is an RNA- 
binding domain (Siomi et al. 1994; Valverde et al. 2008). The 
adaptive substitutions in the Tudor domain suggested their 
possible roles in modulating protein-protein interaction in 
the piRNA pathway. In contrast, LRTs were not significant 
and no sites have been detected in Ago2, Ago3a, Ago3b, 
and Ago4, whereas only two positively selected sites have 
been detected in Ago1. 

3D structure homology modeling further showed that 
some positively selected sites in PiwiM -PAZ and Piwil2-PAZ 
were distributed on pocket of Piwi-PAZ for binding piRNA 
3 r -end, which may facilitate piRNA recognition (fig. 44 and 
B). Some positively selected sites in Tdrdl -third Tudor domain 
and Tdrkh-Tudor domain were localized in the pocket of 
Tudor for binding Piwi-sDMA (fig. 4C and D) to modulate 



Piwi-Tudor protein interaction in the piRNA pathway. These 
analysis highlighted importance of amino acid changes in PAZ 
domain and Tudor domain, which are involved in protein- 
protein/RNA interaction. Overall, these results suggested 
that evolution of functional domains in piRNA pathway 
genes is likely to be adaptive to their functions in piRNA path- 
way in the teleost fish species. 

Higher Substitution Rate in piRNA Pathway Genes in the 
Swamp Eel 

The longest branch in the swamp eel (hermaphroditic fish) 
PiwiM as observed in the teleost lineages (fig. 1) indicated a 
higher substitution in swamp eel than in other fish species. To 
test possible accelerated evolution in this species, we first 
cloned five core piRNA pathway genes (PiwiH, Piwil2, Tdrdl , 
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Fig. 6. — Association between rapid evolution of piRNA pathway and TEs diversity in the teleost fish species. (A) Comparison of TE superfamilies with 
copy number > 1,000 between teleost fish species and mammals. Blue columns, retrotransposon superfamily numbers (Class I); red columns, DNA trans- 
poson superfamily numbers (Class II). The data of human and mouse are from a previous study results (Lander et al. 2001; Mouse Genome Sequencing 
Consortium et al. 2002). As many as 12-25 superfamilies of retrotransposons have been observed in the teleost fish species, whereas only 10-12 super- 
families are detected in the genomes of both human and mouse. For Class II, 7-18 superfamilies are detected in the fish species, whereas there are 2-3 

(continued) 
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Tdrkh, and PrmtS) and five Ago genes (/Ago 7, Ago2, Ago3a, 
Ago3b, and Ago4) from swamp eel using degenerate PCR 
and RACE technique. RT-PCR showed that piRNA pathway 
genes were mainly expressed in three kinds of gonads (ovary, 
ovotestis, and testis), whereas Ago genes were expressed in all 
adult tissues detected, indicating a role of piRNA pathway 
genes in the germline (fig. 5A). We then examined evolution- 
ary rate of piRNA pathway genes in the species. The branch 
model analysis (one-ratio model vs. two-ratio model) showed 
that five piRNA pathway genes (PiwM, Piwil2, Tdrdl, Ddx4, 
and PrmtS) in the swamp eel evolved more rapidly than that in 
other fish species significantly (fig. SB). These results sug- 
gested that piRNA pathway genes in the swamp eel evolved 
at an elevated rate compared with the other gonochoristic 
fishes. 

Diversity of TEs and TE-Derived piRNAs in Teleosts 

To connect rapid evolution of piRNA pathway genes with TEs 
diversity in the teleost fish species, using newly released whole 
genome sequence data, together with new genome (swamp 
eel), TEs in six teleost fish species were analyzed. There were 
12-25 superfamilies of retrotransposons (Class I TEs) detected 
in different fish species (fig. 6A), whereas only 10-12 super- 
families in human and mouse genome. For Class II, 7-18 
superfamilies were detected in the fish species, whereas 
there were 2-3 superfamilies in both human and mouse 
(Lander et al. 2001; Mouse Genome Sequencing 
Consortium et al. 2002). Thus, the teleost fish genomes 
have much higher TE diversity compared with mammalian 
genomes, suggesting a correlation between TEs diversity 
and the rapid evolution of piRNA pathway. Further analysis 
of loci distribution of piRNAs in zebrafish genome showed 
that TE-derived piRNAs were corresponding to 16 superfami- 
lies of retrotransposons and 13 superfamilies of DNA 
transposons (fig. 6B). However, most TE-derived piRNAs 
were originated from about three superfamilies of 
retrotransposons in mouse (Aravin et al. 2007). In addition, 
DNA trsansposons (Class II) are inactive in mammals (Huang 
et al. 2012). These data suggested an association between 
rapid evolution of piRNA pathway and TEs diversity in the 
teleost fish species. 



Discussion 

Advantageous mutation accumulation through positive selec- 
tion may result in adaptive evolution. Adaptive evolution often 
drives rapid evolution (Swanson etal. 2001). Genes involved in 
some biological processes often evolve rapidly. For example, 
immune-related genes are generally thought as a class of rap- 
idly evolving genes, because they are often involved in host- 
parasite revolutionary arms race. In fact, many reproductive 
genes, such as those involved in sex determination (Whitfield 
et al. 1993; Zhang 2004) and gamete recognition (Swanson 
et al. 2001), also evolve rapidly as a result of adaptive evolu- 
tion. In this report, we present that adaptive evolution may 
have shaped piRNA pathway in the teleost fish lineage, and 
the rapid evolution of piRNA pathway is consistent with the 
diversity of TEs, which support a correlation between piRNA 
pathway evolution and transposon diversity in the teleost fish. 
Based on these studies, a model of piRNA pathway in silencing 
diverse TEs in teleost fish genome is proposed (fig. 60- Diverse 
piRNAs transcribed from corresponding TE loci are loaded 
onto Piwi proteins, which have undergone an adaptive evolu- 
tion in the lineage of the teleost fish. Besides Piwi proteins, 
other factors, such as Tdrdl and Tdrkh, are also involved in 
piRNA pathway. Finally, piRNA pathway represses TE activa- 
tion through DNA methylation to protect genome in germline 
cells. 

Studies in Drosophila have indicated that RNAi genes evolve 
rapidly (Obbard et al. 2006; Kolaczkowski et al. 2011). 
Because RNAi processing is associated with TEs silencing, the 
synergetic effects between RNAi and TEs should be important 
in maintenance of genome stability through adaptive evolu- 
tion between small RNA pathway and TEs. Our observations in 
the teleost fish lineages further suggest that rapid evolution of 
fish piRNA pathway is likely to be involved in the adaption to 
transposon diversity. Rapid evolution of piRNA pathway in fish 
may silence TEs more efficiently, which is consistent with ob- 
servations in the Drosophila (Fablet et al. 2014). We find that 
positive selection has occurred in the teleost lineages with 
strongest at the interface between silencing machinery and 
target piRNAs. These results extend our understanding of 
piRNA pathway adaptation to ensure genome stability. 

piRNA pathway genes play crucial roles in the transposon 
repression and the germline functions. Rapid evolution of 



Fig. 6. — Continued 

superfamilies in human and mouse, which are inactive in mammals. The red dots indicate the branches with rapid evolution. (B) Diverse TE loci of zebrafish 
piRNAs. Nonredundant piRNA sequences of zebrafish from Ketting's study (Houwing et al. 2007) together with our data (Zhou et al. 2010) were used to 
search against TE superfamilies using Repeat-Masking with default parameters in Repbase. The number indicates percentage for each superfamily with 
respect to total TE-associated piRNAs. Different colors in pie diagram represent different kinds of the TE superfamily. TE superfamilies with percentages less 
than 1 % are included in other TEs. DIRS, named for the founding member from Dictyostelium discoideum; ERV, endogenous retrovirus; LTR, long terminal 
repeat; L1 , long interspersed element 1 ; L2, long interspersed element 2; non-LTR, nonlong terminal repeat; SINE, short interspersed elements. (Q A model of 
adaptive evolution of piRNA pathway to silence diverse TEs in fish genome. Diverse piRNAs transcribed from corresponding TE loci are loaded onto Piwi 
proteins, which have undergone an adaptive evolution in the lineage of the teleost fish. Besides Piwi proteins, other factors, such as Tdrdl and Tdrkh, are also 
involved in piRNA pathway. Finally, piRNA pathway represses TE activation through DNA methylation to protect genome in germline cells. 
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piRNA pathway genes in the teleost fish lineages is consistent 
with observations in the Drosophila species (Simkin et al. 
2013), supporting an evolutionary "arms race" between TEs 
and the piRNA pathway. Moreover, key functional domains of 
Piwi proteins we detected may serve as hotspots of adaptive 
evolution. The adaptation at MID and PAZ domains of PiwiM 
and Piwil2 provides access to their diverse target piRNA mol- 
ecules, which is consistent with much more diversified TEs 
in the teleost genomes. As an epigenetic regulation, the 
piRNA/TE system is important for transposon silencing to 
ensure integrity of animal genomes. In addition, other epige- 
netic modes are also implicated in TE regulation, such as chro- 
matin structure and DNA methylation (Puszyk et al. 2013). All 
these epigenetic regulations play crucial roles in genome 
defense. 

A distinct biological feature that may shape susceptibility 
to TEs is fertilization mode (Huang et al. 2012). The mam- 
mals reproduce by internal fertilization, an adaptation to 
live on land rather than in water as teleosts do. In contrast, 
most teleosts usually reproduce by external fertilization. As 
both sperm and eggs are mobile in a water environment, 
the fertilization in water is subjected to pressure from harsh 
aquatic environment, climate change, and water pollution. 
In order to ensure sufficient fertilized eggs for species sur- 
vival, external fertilization animals definitely have to shed 
much more sperms and eggs into water than internal fer- 
tilization animals do. Thus, gametogenesis is very active so 
as to maintain sustained reproductive capacity, giving rise 
to higher risk of transposon jumping in genomes at the 
same time. The piRNA pathway genes, as guardians of 
the genome playing crucial role in silencing transposons, 
are supposed to be subjected to strong and persistent se- 
lective pressure in the teleosts. A typical example of the 
rapid evolution in the piRNA pathway is observed in the 
hermaphroditic fish, swamp eel; piRNA pathway genes 
evolved at an elevated rate compared with the other gono- 
choristic fish species. In response to selective pressures 
from diversified TEs and reproductive mode of external fer- 
tilization, piRNA pathway genes in the teleosts are sup- 
posed to evolve faster than their orthologous genes in the 
mammals. 

Conclusions 

The piRNA pathway is responsible for germline specification, 
gametogenesis, transposon silencing, and genome integrity. 
Mobile elements can disrupt genome and its functions. piRNA 
pathway evolution and its adaptation to transposon diversity 
in the teleost fish remain unknown. In this report, we present 
that adaptive evolution is likely to shape piRNA pathway in the 
teleost fish lineages, and the rapid evolution of piRNA path- 
way is associated with the diversity of TEs, which suggested a 
link of adaptation evolution of piRNA pathway to transposon 
diversity in the teleost fish. 
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